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Abstract 
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criterion  was  a  reasonable  indicator  of  the  onset  of  crack  initiation,  the  crack  growth  direction 
was  not  well  predicted,  primarily  due  to  the  development  of  damage  near  the  crack  tip.  In  an  ef¬ 
fort  to  overcome  this  limitation,  a  cohesive  zone  damage  model  was  introduced  to  represent  the 
crack.  This  model  was  then  embedded  in  a  boundary  element  method  so  that  arbitrary  crack 
growth  under  mixed  mode  loading  can  be  simulated.  The  details  of  both  the  cohesive  zone  repre¬ 
sentation  and  the  boundary  element  implementation  are  provided.  Comparison  of  the  predictions 
from  the  stress  intensity  factor  based  fracture  criterion  and  cohesive  zone  model  based  simula¬ 
tion  is  presented. 
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1.  Introduction 


The  mixed-mode  problem  considered  here  involves  determination  of  the  critical  conditions  at 
which  a  crack  subjected  to  a  combination  of  loading  modes  I  and  II  will  initiate.  Another  impor¬ 
tant  consideration  is  the  crack  trajectory  upon  initiation.  Since  the  early  work  of  Erdogan  and  Sih 
(1962),  this  problem  has  received  considerable  attention  (for  example,  Williams  and  Ewing, 
(1972);  Maiti  and  Smith  (1983),  Ueda  et  al.,  (1983))  and  a  number  of  criteria  have  been  pro¬ 
posed  to  predict  crack  initiation  behavior  under  mixed-mode  loading.  In  the  two-dimensional 
case,  the  Cartesian  components  of  the  crack  tip  stresses  can  be  expressed  (in  terms  of  the  polar 
coordinates  as) 

^aP(r,e)  =  -^=fap{B)+^=gap{0)-ir^o,a,p  =  1,2,  (1) 

V2^r  V2;zr 


where  faP(0 )  and  gap{Q)  are  characteristic  functions  and  Kj  and  Kji  are  the  mode  I  and  mode 
II  stress  intensity  factors  respectively.  The  asymmetry  of  the  loading  is  generally  described  in 
terms  of  the  mixity  parameter  ju  which  is  defined  as 


(2) 


Thus  the  mixed-mode  fracture  problem  involves  determining  critical  stress  intensity  factor  pairs 
(kCj,Kj,  )  for  different  values  of  p,  at  which  the  crack  will  initiate.  Furthermore  it  is  required  to 

determine  the  initiation  angle  y  as  a  function  of  p,  as  well  as  the  subsequent  crack  propagation 
path.  Many  different  criteria  have  been  proposed  to  predict  the  critical  condition  for  crack  initia¬ 
tion  and  the  crack  trajectory  upon  initiation;  the  maximum  tangential  stress  criterion  (MTS)  de¬ 
scribed  is  probably  the  most  widely  used.  This  criterion  was  proposed  by  Erdogan  and  Sih 
(1962)  and  is  stated  as  follows:  Crack  extension  starts  at  the  crack  tip  in  a  radial  direction.  This 
extension  is  in  that  radial  direction  perpendicular  to  the  direction  of  the  greatest  tension.  Crack 
extension  begins  when  this  tension  reaches  a  certain  critical  value  at  a  certain  distance  from  the 
crack  tip.  Erdogan  and  Sih  (1962)  used  a  plate  containing  a  central  crack  of  length  2 a  subjected 
to  uniaxial  far- field  loading  a,  the  crack  being  oriented  at  angle  ft  to  the  loading  direction.  They 
obtained  a  relation  for  the  initiation  angle  y  as  a  function  of  p  by  equating  the  partial  derivative 
da ee  /  dO  to  zero  (equivalently  the  shear  stress  ar6  can  be  equated  to  zero).  Both  these  condi¬ 
tions  yield  the  relation 

sin  y  +  (3  cos  y  -  l)p  =  0.  (3) 

From  Eq.(3)  it  can  be  shown  that  for  pure  mode  II  conditions  (ju  =  0)  the  initiation  angle  is  pre¬ 
dicted  to  be  70.5°. 


In  the  present  work,  this  problem  of  mixed-mode  fracture  has  been  examined  for  the  solid 
propellant  materials  which  are  particulate  composites.  In  characterizing  the  fracture  of  solid  pro¬ 
pellants,  theories  of  linear  elastic  fracture  mechanics  (LEFM)  and  linear  and  nonlinear  vis- 
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coelastic  fracture  mechanics  (VEFM)  are  usually  applied.  The  fundamental  constitutive  assump¬ 
tions  of  these  theories  dictate  that  the  theories  will  have  a  certain  range  of  validity  depending  on 
the  nature  of  the  loading,  time  scale  of  application  etc.  For  example,  if  the  loads  are  applied  rap¬ 
idly,  it  might  be  appropriate  to  use  LEFM;  on  the  other  hand,  under  conditions  of  storage,  over 
long  periods  of  time,  viscoelastic  deformation  dominates  and  the  appropriate  theory  to  use  in 
determining  the  reliability  of  the  solid  propellant  would  be  the  VEFM.  While  these  theories  are 
firmly  grounded  in  the  classical  principles  of  mechanics,  experimental  results  obtained  on  solid 
propellant  materials  do  not  appear  to  correlate  well  with  the  predictions  of  the  theory  (Francis  et 
al.,  1980).  The  reasons  for  this  disagreement  could  arise  potentially  from  errors  in  experimental 
measurements  and  in  the  theoretical  interpretation.  The  measurements  typically  obtained  in  the 
experiments  are  of  the  loads  and  displacements  at  the  boundaries  of  the  specimen.  These  are  then 
interpreted  in  terms  of  the  crack  tip  parameters  using  LEFM  or  VEFM  assuming,  of  course,  that 
either  of  these  theories  is  appropriate.  However,  there  are  a  number  of  problems  in  such  an  inter¬ 
pretation.  Near  the  crack  tip,  in  a  small  zone  surrounding  the  crack  tip,  there  exists  a  zone  of 
large  deformations  and  damage,  which  cannot  be  modeled  as  an  elastic  or  viscoelastic  medium. 

We  report  here,  the  results  of  an  investigation  into  mixed-mode  fracture.  In  Section  2,  the 
experimental  methodology  and  results  on  the  experimental  characterization  of  mixed-mode 
fracture  are  presented.  The  problems  associated  with  the  prediction  of  mixed  mode  fracture  are 
described  briefly.  In  an  effort  to  provide  a  predictive  capability  for  describing  mixed-mode  frac¬ 
ture,  the  crack  tip  process  zone  is  modeled  using  a  cohesive  zone  description.  The  idealization  of 
the  crack  tip  processes  in  terms  of  the  cohesive  zone  model  is  described  in  Section  3.  The  cohe¬ 
sive  zone  model  may  be  incorporated  into  any  numerical  scheme  for  generating  crack  growth 
computations,  but  the  boundary  integral  formulation  is  the  most  advantageous  since  it  does  not 
require  extensive  remeshing;  the  BEM  formulation  is  described  in  Section  4.  Two  example 
problems  describing  the  application  of  the  method  to  crack  growth  predictions  are  described  in 
Section  5  in  order  to  demonstrate  the  potential  of  the  method.  Finally,  in  Section  6,  the  results  of 
the  cohesive  zone  model  for  the  solid  propellant  mixed-mode  fracture  are  presented. 


2.  Experiments  on  the  mixed-mode  fracture  of  solid  propellants 

We  examined  the  use  of  a  compact-tension-shear  specimen  to  investigate  the  mixed  mode  frac¬ 
ture  criteria  for  solid  propellants.  This  specimen,  shown  in  Fig.  1  along  with  the  loading  grips, 
was  introduced  by  Buchholz  et  al  (1987)  and  has  been  successfully  in  our  laboratory  in  the  past 
in  examining  mixed  mode  fracture  in  polymers  (Mahajan  and  Ravi-Chandar,  (1989))  and  poly¬ 
mer  interface  cracks  (Foltyn  and  Ravi-Chandar,  (1990)).  The  stress  intensity  factors  Kj  and  Kjj 

for  this  geometry  are  given  as  (see  Buchholz  et  al  (1987)): 
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where  a  is  the  angle  between  the  normal  to  the  crack  and  the  load  line,  P  is  the  applied  load,  a  is 
the  crack  length,  w  is  the  specimen  width  and  t  is  the  specimen  thickness.  In  the  tests  performed 
for  this  report,  these  dimensions  were:  a  =  2in;w  =  4  in  and  t  =  0.2  in.  The  crack  was  introduced 
into  the  specimen  by  cutting  a  slit  to  half  the  width  of  the  specimen.  The  specimen  and  loading 
attachments  were  placed  on  a  specially  designed  test  frame  capable  of  cross  head  rates  of  up  to 
100  in/min;  the  crack  tip  region  was  captured  using  a  videocamera  in  order  to  determine  the  on¬ 
set  of  crack  growth.  The  load  and  load-point-displacement  were  recorded  from  the  test  machine. 

The  experiments  on  the  solid  propellant  specimens  provided  by  the  Phillips  Laboratory  were 
performed  under  constant  rates  of  cross-head  extension;  in  the  range  from  0.01  in/min  and  10 
in/min;  the  angle  a  was  varied  in  the  range  from  0  to  90.  Three  different  temperatures  (-65  F, 
+70  F  or  165  F)  were  considered.  The  details  of  the  experiments  and  their  results  are  reported  in 
the  Final  Report  to  Raytheon  STX  Inc,  who  provided  the  contract  from  Phillips  Laboratory  for 
the  performance  of  the  tests.  Figs.  2,  3  and  4  summarize  the  results  corresponding  to  +70  F.  The 
onset  of  mixed-mode  crack  initiation  is  well  characterized  by  the  MTS  theory  as  seen  in  Fig.  2. 
There  is  a  significant  rate-dependence  for  the  mode  I  stress  intensity  factor  as  shown  in  Fig.  3; 
this  implies  that  a  viscoelastic  model  is  necessary  for  the  characterization  of  the  fracture  of  the 
solid  propellant.  Finally,  in  Fig.  4,  the  variation  of  the  kink  angle  with  the  loading  mixity  is 
shown;  clearly  the  MTS  theory  does  not  predict  the  kink  direction,  especially  at  large  values  of 
mode  II  loading.  The  large  damage  that  is  introduced  near  the  crack  tip  in  the  process  zone  is  be¬ 
lieved  to  be  the  main  reason  for  this  discrepancy.  In  an  effort  to  overcome  this,  a  cohesive  zone 
model  is  proposed  to  represent  the  damaged  zone  near  the  crack  tip  as  described  in  the  next  sec¬ 
tion. 


3.  Cohesive  zone  model  for  a  crack  in  elastostatics 

In  this  section,  we  consider  the  idealized  model  of  a  line  fracture  process  zone  near  the  crack  tip. 
This  model  of  the  fracture  process  zone  is  motivated  by  the  fact  that  in  some  materials  the  crack 
surfaces  are  usually  not  separated  completely  behind  the  (fictitious)  crack  tip.  There  exists  a 
relatively  long  extension  of  the  crack  -  variously  called  the  wake  zone,  the  bridging  zone,  or  co¬ 
hesive  zone  -  where  tractions  can  be  transferred  across  the  crack  line.  In  solid  propellants,  this  is 
the  damage  zone  in  the  crack  tip  region.  The  key  assumption  in  this  model  is  that  material  sof¬ 
tening  or  damage  accumulation  beyond  the  peak  load  is  localized  in  a  narrow  layer  behind  a  fic¬ 
titious  crack  tip,  whose  volume  is  negligible  and  whose  action  is  replaceable  by  cohesive  forces. 
Typically,  two  types  of  constitutive  laws  are  used  in  the  literature  for  cohesive  materials:  one  is 
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characterized  by  a  traction-displacement  relationship,  the  other  one  by  a  material  constitutive  law 
defined  in  terms  of  stress  and  strain  accompanied  with  a  thickening  law  of  the  layer.  In  the  latter, 
thickening  of  a  cohesive  layer  is  decoupled  which  may  be  worthy  of  studying  separately  in  vari¬ 
ous  materials  considering  rate-dependence  and  under  dynamic  loading.  However,  for  quasistatic 
loading  -  the  case  considered  predominantly  in  the  literature  -  a  law  describing  the  traction- 
displacement  relation  is  sufficient  for  modeling  of  the  cohesive  zone.  We  present  below  one  such 
law  for  a  cohesive  crack  represented  by  a  line  spring  connecting  two  coincident  crack  points  - 
this  seems  to  be  the  most  appropriate  model  for  crazing-dominant  brittle  polymers  and  fiber- 
reinforced  composites. 

Consider  a  representation  of  the  development  of  a  crack  shown  in  Fig.  5a.  We  suggest  that 
two  points,  x+  and  x" ,  originally  coincident  on  opposite  sides  of  a  line,  separate  into  two  dis¬ 
tinct  points,  connected  by  the  cohesive  zone  material;  continued  straining  increases  the  separa¬ 
tion  between  these  two  points  and  eventually  leads  to  cracking.  The  kinematics  of  this  separation 
process  are  assumed  to  be  described  completely  by  the  crack  face  separation,  w.  Introducing  the 
local  normal  and  tangent  directions  along  the  cohesive  zone,  w  can  be  resolved  into  the  normal 
separation  distance  (or  the  cohesive  zone  opening  displacement)  component,  wn  =  u~  -  u*  and 

the  tangential  separation  distance  (or  the  sliding  displacement)  component  wt  =u~  -u+ .  In  or¬ 
der  to  prevent  the  inter-penetration  of  the  cohesive  zone,  wn  >  0 ;  equality  holds  only  in  the  case 

that  there  is  contact  between  the  top  and  bottom  surfaces  of  the  cohesive  zone  and  in  this  case, 
we  must  also  provide  a  description  of  the  frictional  resistance  on  the  surface  as  well.  We  con¬ 
sider  first  the  case  of  a  locally  opening  mode  crack  with  wn  >  0  and  describe  the  force- 

separation  law  for  the  cohesive  zone.  The  cohesive  material  may  be  modeled  by  a  simple  line 
spring  that  behaves  according  to  the  following: 


p  =  k(wrf)w, 

(5) 

or  in  component  form 

Pn  =k(wd)wn 

Pr  =k(wd)wT  ’ 

(6) 

where  p  is  the  traction  vector  with  normal  and  tangential  components,  pn  and  pT ,  respectively. 

wj  is  the  maximum  separation  distance  between  two  originally  coincident  points  on  the  crack 
over  the  entire  loading  history,  and  is  used  as  a  damage  parameter.  The  stiffness  of  the  cohesive 
zone  material  is  denoted  by  k{wj)  and  is  assumed  to  depend  on  the  current  state  of  damage.  Note 
that  k(wd)  is  a  decreasing  function  of  wj,  indicating  the  softening  behavior  of  the  material;  this 
imposes  an  irreversibility  of  the  damage  process  under  unloading.  Corresponding  to  each  Wd, 
there  exists  a  traction,  pd,  and  a  damage  locus  derived  from  Eqs.  (5): 

Pd  =  k(wd)wd  (7) 

The  constitutive  law  described  above  is  illustrated  schematically  in  Fig.  5b.  Note  that  the  de¬ 
scription  of  the  cohesive  zone  material  through  the  damage  parameter  and  stiffness  allows  for 
irreversibility  of  damage.  Upon  unloading,  the  points  on  the  cohesive  zone  unload  linearly  with  a 
stiffness  k(wd),  whereas  in  most  cohesive  zone  models,  unloading  occurs  along  the  damage  locus. 
Including  the  unloading  effect  has  expanded  the  capability  of  the  cohesive  crack  model  in  simu¬ 
lating  realistic  fracture  behavior  significantly.  There  are  two  critical  states  along  the  damage  lo- 
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cus.  The  first  one,  at  wj  =  0  and  pd  =  py,  represents  the  maximum  traction  that  can  be  sustained 
by  the  material  before  the  cohesive  zone  begins  to  develop;  beyond  this  critical  level,  separation 
processes  begin  and  wd  increases.  The  point  on  the  specimen  that  is  at  this  state  is  usually  called 
the  fictitious  crack  tip  or  the  cohesive  zone  tip.  The  second  critical  point  on  the  damage  locus 
occurs  at  wd  =  w/  and  pd  =  0;  this  point  represents  the  maximum  displacement  jump  across  the 
cohesive  zone  that  can  be  sustained  before  cracking;  beyond  this  level,  the  traction  goes  to  zero 
and  the  two  initially  coincident  points  are  now  completely  separated.  The  point  on  the  specimen 
that  is  at  this  state  is  usually  called  the  physical  crack  tip.  Hence,  Eqs.  (5)  and  (6)  together,  define 
the  complete  process  of  separation  of  a  material  point  into  a  crack,  as  long  as  wn  >  0 .  A  consti¬ 
tutive  law  in  stress-strain  for  a  cohesive  material  can  be  described  in  a  similar  manner  (see  Yang 
and  Ravi-Chandar,  1996).  The  transition  of  that  law  to  this  law  is  apparent. 

We  now  turn  to  the  case  when  wn  =  0 .  Under  arbitrary  loading,  contact  of  two  crack  surfaces 

could  occur.  In  such  cases,  as  shown  in  Fig.  6a,  we  must  provide  an  appropriate  description  for 
the  development  of  friction  along  the  cohesive  zone,  in  addition  to  the  line-spring  model  de¬ 
scribed  above.  We  note  that  this  is  a  generic  development  and  not  particularly  applicable  to  the 
case  of  the  solid  propellant;  for  the  propellant,  it  might  be  appropriate  to  set  the  friction  coeffi¬ 
cient  to  zero.  We  model  the  frictional  interaction  simply  through  a  Coulomb  type  law.  Thus,  the 
tangential  interaction  of  the  crack  surfaces  in  contact  is  given  by 

Pt  =  Kwd  K  -  sgn(wr  )fpn ,  (8) 

where /is  frictional  coefficient.  The  equation  for  the  determination  of  the  normal  component  of 
the  traction  pn  is  obtained  by  enforcing  the  contact  condition  wn  =  0.  In  Eq.  (8),  the  first  term 

represents  the  tangential  traction  contribution  due  to  the  line  spring  described  above  (nontrivially 
if  it  is  not  broken  completely)  and  the  second  term  represents  the  Coulomb  friction  component. 
In  order  to  attain  smooth  transition  of  frictional  force  near  wT  =  0,  the  frictional  coefficient / may 
be  assumed  to  be 


l/o> 


K  p  w*; 

otherwise, 


(9) 


where /0  and  are  both  non-negative  material  constants.  Function /is  plotted  in  Fig.  6b.  This 
completes  the  description  of  the  cohesive  zone  material  behavior. 

There  still  remains  the  issue  of  deciding  on  the  appropriate  incorporation  of  this  model  into  an 
elastostatic  crack  problem.  One  major  question  that  arises  is  the  following:  What  is  the  criterion 
that  can  be  used  to  grow  the  cohesive  zone  from  a  stress  contentrator?  The  incorporation  of  the 
cohesive  zone  model  eliminates  energy  release  criterion  from  consideration  since  the  energy  re¬ 
lease  rate  will  always  be  the  same  regardless  of  the  direction  of  crack  extension  -  it  is  simply  the 
area  under  the  damage  locus,  assuming  a  fully  developed  cohesive  zone.  The  most  plausible  cri¬ 
teria,  particularly  for  the  line-spring  nature  of  the  cohesive  zone  mode,  are  stress  based  criteria 
such  as  the  maximum  principal  stress  criterion  or  the  maximum  tangential  stress  criterion  at  the 
fictitious  crack  tip.  Assuming  that  one  of  these  criteria  would  be  appropriate,  a  second  major 
question  arises:  What  should  be  the  step  size  in  extending  the  crack  along  this  direction?  In  finite 
element  formulations,  such  as  those  of  Xu  and  Needleman,  (1994),  and  Ortiz  (1996),  cohesive 
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zones  are  forced  to  develop  along  element  boundaries  and  the  extension  is  over  the  side  of  one 
entire  element.  The  approach  in  these  models  is  to  make  the  element  size  so  small  that  the  overall 
crack  growth  behavior  is  captured  adequately.  In  other  words,  the  macroscopic  crack  path  is  sug¬ 
gested  to  be  independent  of  the  length  scale  of  the  discretization  when  the  latter  is  sufficiently 
small  compared  to  characteristic  structural  length  scale.  In  the  present  work,  we  choose  the 
maximum  principal  stress  criterion  for  determining  the  crack  increment.  If  the  maximum  princi¬ 
pal  stress  at  a  fictitious  crack  tip  reaches  the  critical  value,  py,  this  tip  is  ready  to  run  under  fur¬ 
ther  loading.  The  direction  in  which  the  tip  advances  is  perpendicular  to  the  direction  of  the 
maximum  principal  stress  at  that  point,  and  the  extension  is  such  that  the  maximum  principal 
stress  at  the  new  tip  position  is  kept  at  the  critical  value  py,  during  continued  loading.  However, 
note  that,  structural  instability  may  occur  while  a  crack  is  advancing;  i.e.,  extension  of  a  crack 
may  enhance  the  stress  state  at  a  crack  tip  rather  than  release  it.  In  this  case,  the  crack  may  run 
fast  and  inertia  effects  may  have  to  be  included. 


4.  Boundary  integral  formulation  of  the  elastostatic  crack  problem  with  a  cohesive  zone  at 
the  crack  tip 

Boundary  element  method  (BEM)  has  received  much  attention  recently,  especially  in  application 
to  fracture  mechanics  (FM)  (see  Cruse,  1996  and  Aliabadi,  1997  for  reviews).  The  method  is  at¬ 
tractive  because  it  involves  discretization  of  the  boundary  alone;  the  dimensionality  of  the  stiff¬ 
ness  matrix  formed  in  BEM  is  then  reduced  by  one  in  comparison  to  a  domain  method,  such  as 
finite  element  method  (FEM),  although  the  stiffness  matrix  is  full  and  asymmetric  in  general.  A 
particularly  attractive  advantage  of  BEM  in  application  to  crack  problems  over  FEM  is  that  do¬ 
main  remeshing  is  not  necessary  when  a  crack  grows;  only  one  more  element  of  a  crack  is  added 
with  all  the  already  existing  elements  untouched.  However,  as  a  crack  is  modeled  mathematically 
with  the  two  crack  surfaces  being  coincident,  the  classic  boundary  integral  equation  (BIE)  can 
not  be  applied  directly  or  the  resulting  stiffness  matrix  formed  will  be  ill-conditioned.  A  great 
deal  of  effort  has  been  expended  in  dealing  with  this  difficulty;  many  methods  -  such  as  the  crack 
Green’s  function  method  (Snyder  and  Cruse,  1975),  the  displacement  discontinuity  method 
(Crouch,  1976,  and  Wen,  1996),  the  subdomains  method  (Blandford,  Ingraffea  and  Liggett, 
1981),  the  dual  boundary  element  method  (DBEM)  (Hong  and  Chen,  1988,  Portela,  et  al.,  1992, 
and  Chen  and  Chen,  1995),  the  single-domain  traction  boundary  element  method  (Young,  1996), 
and  some  hybrid  methods  optimizing  the  advantages  of  some  of  the  above  mentioned  methods 
(Ameen  and  Raghuprasad,  1994)  -  have  been  proposed.  It  is  not  our  intent  to  review  these  meth¬ 
ods  here;  one  may  see  Chen  and  Chen  (1995)  for  a  discussion  of  some  of  these  methods.  Among 
the  methods  mentioned  above  the  method  of  subdomains  eliminates  the  discrepancy  noted  above 
by  cutting  a  cracked  structure  into  pieces  of  simpler  topology  such  that  BIE  can  be  applied  prop¬ 
erly  to  each  of  domain  separately.  However,  this  formulation  increases  the  computational  effort 
as  a  result  of  the  additional  artificial  boundaries  that  appear  from  the  cutting  process;  moreover, 
when  a  crack  advances,  remeshing  is  still  needed  in  general.  Nevertheless,  the  idea  of  the  method 
of  subdomains  is  valuable  and  is  adopted  in  the  present  paper  to  derive  a  mathematically  rigor¬ 
ous  and  simple  formulation  of  the  single-domain  dual-boundary-integral  equations  (SDDBIEs) 
of  a  cracked  structure.  The  classic  dual  integral  equations  for  a  simple  structure  is,  of  course, 
well  established  in  the  literature.  These  SDDBIEs  were  given  by  Young  (1996),  and  applied  to 
solve  a  traction-free  crack  problem  using  continuous  elements. 
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4.1.  Summary  of  the  dual  integral  equations  for  a  simple  structure 

Consider  a  homogeneous,  isotropic,  linearly  elastic  domain  Q  with  piecewise  smooth  boundary 
T  as  shown  in  Fig.  7.  A  Cartesian  (Lagrangian)  coordinate  system1  is  used  along  with  standard 
indicial  notation.  The  displacement  components  ut  (X),  at  a  point  X,  can  be  represented  in  the 
following  form: 

c0(X)uj(X)  =  J  { « ‘  (X,  x)P]  ( x)  -  />  *  (X,  x,  n)w;  (x)  }</T(  x)  +  |M*(X,x)6;(x)(iQ(x),  (10) 

r  v  n 


where  P]  are  traction  components,  bj  are  body  force  components.  «*  and  Pij  are  the  fundamen¬ 
tal  solutions  representing  the  displacements  and  tractions  respectively  in  the  f'  direction  at  a 
field  point  x  due  to  a  unit  force  acting  in  the  ilh  direction  at  a  source  point  X.  Note  that  P*  are 

taken  along  the  outward  normal  vector,  n,  of  T  at  x  and  that  with  the  inclusion  of  n  in  P* ,  the 
dependence  of  the  integrals  on  the  normal  is  made  explicit.  Cy  is  a  coefficient  matrix  given  by 


4  XeQ; 
ctj  XeT; 

0  otherwise, 


(ID 


where  8tj  is  Kronecker  delta,  and  Cy  =  Syt. 2  if  the  tangential  surface  at  X  is  smooth;  if  it  is  not  the 

case  one  may  see  Hartmann  (1980)  for  the  closed-form  expressions  of  this  matrix.  When  XeQ, 
Eqs.  (1)  with  Cy  =  5tj  is  called  the  Somigliana  identity.  One  may  easily  prove  that,  if  X  is  outside 

Q  and  T ,  Eqs.  (10)  holds  trivially  with  c,y  =  0.  The  form  of  Eqs.  (10)  most  useful  in  the  bound¬ 
ary  integral  formulation  arises  when  XeT.  Applying  a  limiting  process  as  X  approaches  a 
boundary,  the  Somigliana  identity  leads  to  the  boundary  integral  equation  of  displacements 
(BIED,  or  usually  BIE),  based  on  which  the  classic  boundary  element  technique  is  developed 
(Brebbia,  Telles  and  Wrobel,  1984). 

Differentiating  Eqs.  (10)  with  respect  to  X  (with  X  eQ)  the  strains  over  Q  can  be  calculated. 
If  these  strains  are  substituted  into  Hooke’s  law,  the  integral  equations  of  stresses  at  a  source 
point  X  (with  X  e  Q )  are  obtained.  A  limiting  process,  similar  to  the  one  used  to  obtain  the  BIE, 
can  be  applied  to  the  resulting  integral  equation  as  X  approaches  a  boundary,  leading  to  the 
boundary  integral  equation  of  stresses.  These  integral  equations  are  summarized  below: 

C,(X)<r,(X)  =  J(c/;(X,x)p,(x)-/>;(X,x,nK(x)}<fl'(*)+  (12) 

r  n 

where  U*Jk  and  P‘k  are  linear  combinations  of  derivatives  of  u*  and  p* ,  with  respect  to  X.  As  in 

the  case  of  Plj ,  P*k  are  taken  along  the  outward  normal  vector,  n,  of  T  at  x.  Cy  is  a  coefficient 
matrix  given  by 


1  Standard  index  notation  is  used.  The  range  for  Latin  subscript  indices  is  3  and  the  range  for  Greek  subscript  indices 
is  2.  Summation  over  repeated  subscripts  over  their  range  is  implied  unless  suspended  explicitly.  Superscripts  do  not 
follow  this  range  and  summation  convention;  their  range  would  be  indicated  explicitly  and  summation  is  indicated 
by  the  summation  symbol. 
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X  gQ; 
XeT; 
otherwise, 


(13) 


C 


4 

<4/2 

0 


which  is  identical  to  c,j  provided  that  comer  points  where  tractions  are  not  well  defined  are  ex¬ 
cluded.  If  Eqs.  (12)  are  multiplied  on  both  sides  by  the  outward  normal  vector  at  X  as  XeT, 
then,  the  boundary  integral  equations  of  tractions  (BIET)  are  obtained,  which  may  be  employed 
to  formulate  a  boundary  element  method  in  a  similar  way  as  BIE. 


The  fundamental  solutions  used  in  the  above  formulation  are  due  to  Kelvin  (Love,  1944),  and 
are  the  basic  singular  solutions  to  a  point  load  in  an  infinite  medium.  The  components  u*  and  /?* 
are  given  below  for  two  (plane-strain)  and  three  dimensional  problems: 


"*(X’I)=  v)r  [{3  -  4V)^  +  W]>  (3-D) 


«;(x,x)= 


8^(1- v)L 

j!?‘(X,x,n)  = 


(3-4v)lnQ^+r,r, 


(2-D;  plane-strain) 


-1 

dr 

4an{\-  v)ra 

dn 

x((1_2v/)4  +  Prs,j)-{x-2vfr,<nj  ~rjn) 


(14) 

(15) 

(16) 


where  p  is  shear  modulus,  and  v  is  Poisson  ratio;  a  =  2,  /?  =  3  for  three  dimensional  prob¬ 
lems  and  a  =  1,  /?  =  2  for  two-dimensional  plane-strain  problems.  Moreover,  r  =  r(x,x)  rep¬ 
resents  the  distance  between  the  points  X  and  x,  and  its  derivatives  are  taken  with  respect  to  x. 


Uijk  and  P*k  are  given  by 


+  Pv{rinjrk  +n,rJrJt 


[(i  -  2v)(rj8lk  +  r,dJk  -  r k8y )+  ] , 

j\p^1  ~  2v)r*sij +  v(r./4* + rJ5ik )-  nr,k  ] 
)+  (l  -  2  +  nj8tk  +  n,8jk )-  (1  -  4  v)nk8tj  }, 


(17) 


(18) 


where  a  and  p  are  as  given  above,  and  y  =  5  for  three  dimensional  problems  and  y  —  4  for 
two  dimensional  problems.  The  plane  strain  expressions  are  valid  for  plane  stress  provided  that 
v  is  replaced  by  v  =  v/(l  +  v) .  As  a  source  point  X  and  a  field  point  x  coincide  on  a  boundary, 

the  integral  kernels  of  «*. ,  ptJ ,  U'jk ,  and  P'jk  may  be  singular  and  even  hypersingular.  The  corre¬ 
sponding  integrals  in  Eqs.  (10)  and  (12)  are  taken  in  the  sense  of  Cauchy  principal  value  if  sin¬ 
gular,  and  are  taken  in  the  sense  of  Hadamard  principal  value  if  hypersingular. 


4.2.  Single-domain  dual  integral  equations  for  a  cracked  structure 

A  structure  containing  a  mathematically  sharp  crack  degenerates  the  boundary  integral  formula¬ 
tion  due  to  coincidence  of  the  two  crack  surfaces.  In  that  case,  one  cannot  apply  either  BIED  or 
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BIET  to  produce  a  boundary  element  method  directly  in  general.  This  problem  has  been  ad¬ 
dressed  in  a  number  of  ways:  partitioning  the  domain  into  multi-domains,  using  crack  Green’s 
functions,  the  displacement  discontinuity  technique,  and  so  on.  In  this  section,  we  describe  the 
formulation  of  a  single-domain  dual-boundary-integral  method  using  Kelvin’s  solutions.  We  cut 
a  cracked  elastostatic  structure  into  subdomains  of  simple  topology  so  that  the  dual  integral 
equations  for  a  simple  structure  can  be  applied  to  each  domain  appropriately.  In  a  further  step, 
using  the  continuity  and  equilibrium  conditions  along  artificial  boundaries,  as  well  as  some  prop¬ 
erties  of  the  integral  kernels,  we  eliminate  entirely  the  integrals  involving  any  artificial  bounda¬ 
ries.  In  order  to  illustrate  the  derivation  of  the  single-domain  integral  equations,  we  shall  con¬ 
sider  a  structure  containing  one  crack;  however,  the  same  approach  can  be  applied  to  any  number 
of  cracks  in  the  structure.  Furthermore,  for  simplicity,  we  shall  ignore  the  body  forces;  if  neces¬ 
sary,  they  may  be  added  back  into  the  final  equations  without  any  difficulty. 

As  shown  in  Fig.  8,  a  structure  containing  one  crack  is  cut  into  two  subdomains  with  a  path 
that  passes  through  the  crack.  We  denote  the  two  subdomains  without  a  crack  by  Q,  and  fin 
respectively.  Their  boundaries  are  denoted  by  Tj  and  Tn  respectively.  The  topology  of  the  sub- 
domains  is  simple  such  that  Eqs.  (10)  and  (12)  may  be  applied  to  them.  This  results  in  the  fol¬ 
lowing  sets  of  equations  for  the  displacement  and  traction  components  in  the  two  subdomains. 
For  the  domain  Q, ,  we  have 

4(X)«;.(X)=  /{ M-  (X,  (x)  -  (X,  x,  n1  )»>  (x)}^r(x) ,  (19) 

b 

with 


cl 

-ij 


4 

4 


0 


X  eQ,; 

Xer,; 

otherwise. 


(20) 


and 


C|,(XK(X)  =  J{l/;(X,xK(x)  -  ^(X,x,n‘>;(x)}<r(*) , 
r, 

with 


Cl 


Sik  X  ef2j; 

<4/2  X€Ti; 

0  otherwise 


(21) 


(22) 


For  the  domain  fiu ,  we  also  have 

c“(X)k,"(X)=  J{U;(X,x)^,(x)-A;(x,x.nll)»»(x)j<r(x),  (23) 

r„ v 

with 


9 


=  < 


(24) 


4- 


Sy  X  €Qn; 

c"  Xeru; 

0  otherwise, 


and 


£!(x)4(x)  =  J{c/;(x,i)A"(*)-p;(x,x,n")»;(i)}dr(x), 


(25) 


with 


Cl  = 


^ik 


8ik,  X  e  Qn ; 

4/2,  X€r„; 


(26) 


0, 


otherwise. 


All  of  the  functions  are  defined  for  all  X  and  are  single-valued.  Thus,  we  add  Eqs.  (19)  and  (23), 
and  Eqs.  (21)  and  (25)  algebraically,  and  rearrange  them  in  the  following  forms: 


4  (X)m'  ( x)  +  c  J  (X)w“  ( X)  =  J \u]j  (X,  x)pj  (x)  -  p*  (X,  x,  n)uj  ( x)  |g?T(  x)  + 

rex 

j  { 4 (X, x)[4 ( x)  +  />" ( x)] - [p * (x, x, n 1  )w ‘ ( x)  +  p *(x, X, n 11  )ulJ ( x)  ]  W(  x) , 


(27) 


and 


WoJ (X) + c“  (x)<(x)-  j{u;„ (x,\)Pk(x)  -  p;„ (x,x,n K (x)}dr(x) + 
f  { u'jk (X, x)[/?'  (x)  +  pf  (x)]  -  [i£  (x, X, n1  )u\  (x)  +  p;k  (x,x,nn  )w" (x)]j dT{x) , 


(28) 


where  the  superscript  ex  denotes  an  external  regular  boundary,  the  superscript  c  denotes  a  crack 
boundary,  and  the  superscript  a  denotes  an  artificial  boundary  generated  by  the  cutting  process. 
Moreover,  the  integrals  over  T“  and  T,f  have  been  put  together  with  T“  =  Tj“  +  T,f ;  the  su¬ 
perscripts  I  and  II  have  been  dropped  since  T“  and  Tj“  do  not  overlap.  Equilibrium  of  the 
whole  structure  indicates  that  the  tractions  on  two  coincident  opposite  artificial  boundaries  are 


equal  in  magnitude  and  opposite  in  direction;  i.e., 

p){x)  =  -pf{x),  as  x  e  T°.  (29) 

Also,  due  to  the  requirement  of  continuity  of  the  displacement  components  at  the  coincident 
points  on  the  artificial  boundaries,  it  holds  that 

u,‘(x)  =  u!I(x),  as  x  €  r°.  (30) 

Moreover,  the  integral  kernels  in  Eqs.  (7)  and  (9)  have  the  following  properties: 

4(x>x  ~n) =  ~Po(x>x>n)>  i  e->  />*(X>XV)  =  -4(x,x,nn)  as  x  e  rc+a ,  (31) 

P*k(x,x,-n)  =  -P*k(x,x,n),i.e.,  P^(x,x,nl)=  -P’k(x,x,n11)  asx  €  rc+°.  (32) 
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In  deriving  the  above,  the  reversal  of  the  outward  normal  n  at  the  two  coincident  points  on  the 
crack  and  artificial  boundaries  has  been  applied: 

n'(x)  =  -nn(x)  asx  eP.  (33) 


Using  Eqs.  (15),  (17),  and  (20)  through  (23)  to  Eqs.  (18)  and  (19),  we  obtain 


7,(x)  =  {{M’(x,x)p7(x)-jp*(x,x,n)M7(x)}jr(x) 

rex 

+  J  { « *■  (X,  x)(pj  (x)  +  pj(x))~  p*  (x,  x,  n + )(«;  ( x)  -  (x)|  dr{  x) , 

rc 


with 


/,(x) 


m,(X)  Xefi; 

^  cff(X)uj(X)  XeP; 

\;(x)w;(x)+c-(x)M;(x)  x€r; 

0  otherwise, 


(34) 


(35) 


and 


*4(X)  =  f{u*t(X,x)pt(x)  -  i>;(X,x,n)M,(x)}rfT(x) 

rex 

+  \{ulk (X, \)(p+k  (x)  +  p~k{x)) -  P’k (x, x, n+ \uk (x)  - uk  (x))}</T(x) , 
rc 

with 


S,(X) 


'ov(X) 

<o-,(X)/2 

'(cr*(X)  +  <T-(X))/2 
0 


X  eQ; 
XeP; 
XeTc; 
otherwise, 


(36) 


(37) 


where  the  superscript  +  indicates  one  side  of  a  crack,  and  the  superscript '  indicates  the  corre¬ 
sponding  opposite  side,  instead  of  the  superscripts  I  and  II.  The  positive  side  of  a  crack  may  be 
chosen  arbitrarily,  for  convenience.  Eqs.  (34)  and  (36)  are  so-called  the  single-domain  dual  inte¬ 
gral  equations  of  a  cracked  structure ;  in  these  equations  the  integrals  are  taken  only  along  the 
regular  external  boundary  and  one  side  of  a  crack.  The  artificial  boundaries  that  appeared  due  to 
the  cutting  process  have  been  entirely  eliminated  from  consideration.  If  a  structure  contains  mul¬ 
tiple  cracks,  the  single-domain  dual  integral  equations  are  given  in  the  same  forms  as  Eqs.  (34) 
and  (36). 


Eqs.  (34)  and  (36)  with  a  source  point  on  the  boundary,  i.e.,  X  eT“+c,  are  of  the  most  im¬ 
portance,  in  the  formulation  of  boundary  integral  equations  for  elastostatic  problems.  In  formu¬ 
lating  a  general  boundary  integral  method  using  Eqs.  (36),  it  is  advantageous  to  use  tractions  in¬ 
stead  of  stresses  since  the  number  of  unknowns  can  be  reduced.  If  Eqs.  (36)  are  “multiplied”  on 
both  sides  by  the  outward  normal  at  the  boundary  point  X  as  usual  the  traction  version  of  Eqs. 
(36)  with  X  e  Tex+C  can  be  obtained  as 
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(38) 


4<xK(x)  =  |{t/;(x,x)pi(x)-j»;(x,i,ii)»4(i)}i,>(x)<r(i) 

rer 

+  j { ulk  (x>  *)(p+k  ( x)  +  P,  ( x) )  -  />;  (x, x,  n +  X«;  ( x)  -  w;  ( x) )} n}  ( X)  rfT(  x) . 

rc 

Note  that,  as  X  eTc,  n/X)  is  taken  to  be  the  outward  normal  of  the  positive  crack  side,  i.e., 

hXX)=»;(x). 

Note  that  on  the  crack  surfaces,  the  tractions  are  self-equilibrating  and  hence 

»;(X)(<r;(X)  +  cr:(X))/2  =  A*(X)  and  p,*(x)  +  p-(x)  =  0.  (39) 

Eqs.  (39)  can  be  used  to  simplify  Eqs.  (34),  (36)  and  (38)  when  applying  to  cracks.  In  this  case, 
traction  and  displacement  jump  on  a  crack  may  be  regarded  as  independent  variables  in  a  nu¬ 
merical  formulation,  with  the  least  unknowns. 

4.3.  Iterative  boundary  element  method  of  successive  over-relaxation 

The  behavior  of  the  cohesive  zone  materials  and  the  criterion  for  crack  advance  described  in  the 
last  section  require  that  an  incremental  loading  procedure  be  used  when  solving  boundary  value 
problems.  Moreover,  an  iteration  process  is  indispensable  in  each  loading  step  in  general  due  to 
the  fact  that  the  constitutive  law  of  the  cohesive  zone  material,  given  in  Eqs.  (5)  -  (9),  involves 
irreversible  damaging  process,  and  is  essentially  history-dependent.  Furthermore,  the  stiffness  of 
the  model  line  spring  of  the  cohesive  zone  material,  Jdwj),  is  in  general  a  part  of  the  solution, 
which  may  or  may  not  be  dependent  on  the  current  displacement  discontinuity;  in  other  words, 
there  is  no  simple  relationship  between  the  traction  and  displacement  discontinuity  on  a  crack 
which  may  be  used  to  achieve  a  linear  system  of  equations  to  solve  the  problems  numerically. 
Two  strategies  of  iteration  are  possible  in  general:  in  the  first,  a  linear  system  of  equations  of  the 
discretized  problem  is  formulated  using  the  stiffness  of  the  cohesive  zone  material  obtained  in 
the  previous  iteration  step,  and  is  solved  by  a  typical  solver  either  direct  or  iterative.  The  stiffness 
of  the  cohesive  zone  material  is  then  modified  based  on  this  solution,  and  the  process  is  repeated 
until  the  solution  of  the  desired  accuracy  is  achieved.  In  the  second  strategy,  a  nonlinear  system 
of  equations  of  the  discretized  problem  is  formulated  using  the  stiffness  of  the  cohesive  zone 
material  as  unknowns,  and  is  solved  iteratively.  By  the  first  strategy,  the  iterative  procedure  is 
very  clear  and  the  convergence  of  the  iteration  process  is  expected.  However,  it  is  apparently 
time-consuming  even  with  a  very  efficient  solver  of  a  linear  system  of  equations;  the  linear  sys¬ 
tem  of  equations  is  solved  fully  many  times  until  a  solution  is  achieved.  On  the  other  hand,  the 
total  time  for  a  solution  by  the  second  strategy,  if  a  good  iterative  procedure  is  found,  would  be 
in  the  same  order  as  that  for  one  step  of  the  iteration  by  the  first  strategy  (solving  the  linear  sys¬ 
tem  of  equations  once).  In  the  present  work,  we  adopt  the  second  strategy  and  formulate  an  it¬ 
erative  method  of  successive  over-relaxation  for  the  present  nonlinear  problem  of  a  cohesive 
crack.  The  solution  procedure  of  this  iterative  method  will  be  described  below,  following  discre¬ 
tization  of  the  boundaries  and  the  boundary  integral  equations. 

The  boundaries  of  a  cracked  two-dimensional  structure  are  approximated  by  straight  elements, 
Tel  s,  of  which  each  contains  Nei  nodes  that  are  uniformly  distributed  in  it,  as  shown  in  Fig.  9. 
These  nodes  are  numbered  separately  on  the  external  boundary  and  on  the  crack  locus.  Assume 
that  we  have  Ifx  nodes  on  the  external  boundary  and  If  nodes  on  the  crack  locus.  A  field  quan- 
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tity,  q(x)  can  be  approximated  over  an  element  Td  by  interpolating  the  nodal  values,  q"  in  this 
element  as 


q(x)  =  fy(x)q\ 

«  =  1 


(40) 


where  the  interpolation  function,  ^”(x)  satisfies  the  following  conditions: 


m  =  n; 
m^n, 


and  Yjf{\)  = 


fl  xeTel ; 

jo  xerel. 


(41) 


The  interpolation  functions  may  be  constant,  linear  or  quadratic,  depending  on  the  required  accu¬ 
racy  of  the  representation.  Substitute  Eq.  (40)  into  Eqs.  (34)  and  (38)  with  X  e  Fex+C ,  the  discre¬ 
tized  forms  of  the  dual  boundary  integral  equations  are  obtained.  Thus,  one  obtains 

Na  Na+Nc 

-yx")  =  0,  m  =  (42) 

n=l  n=Nex+ 1 

and 


N"  Na+Nc 

Y,\p:;pnp-Hma;unpy  -W>o,  ^  =  1,2,...^*+^.  (43) 

«=1  n=V"+ 1 

Eqs.  (38)  and  (39)  each  represent  2(A/ex  +  tf  )  equations  that  are  the  discretized  version  of  Eqs 
(25)  and  (29).  In  these  equations,  ,  h"}} ,  ,  and  are  given  explicitly  by: 


gap  =  fMi?(x^xV”(x)<*r(x)> 

T” 

(44) 

Kp  =  J/v(x^x>nV”(xWx)’ 

r” 

(45) 

g:;  =  JtC,(x-.*y’(xK(x-)r(x), 

r" 

(46) 

HTf  =  JC(x^*>“V',(*k>(x')a'M. 

(47) 

r" 


where  T”  is  the  element  where  the  nth  node  is  located.  Eqs.  (44)  through  (47)  may  be  evaluated 
numerically  if  Xm  g  T" ,  and  analytically  if  Xm  eT".  Note  that  the  condition  for  a  self- 
equilibrating  crack,  Eq.  (39),  has  been  used  in  the  above  equations.  If  u  and  p,  are  a  trial  set  of 
displacements  and  tractions,  Eqs.  (42)  and  (43)  would  not  be  satisfied  in  general;  instead,  we 
obtain  the  following  residuals: 


N‘*+Nc 


kk)+ 


n=Nex+ 1 


m  =  l,2,...Nex  +  Nc, 


(48) 


Na+Nc 


-h:k)+ 

n=\  n=Nex+ 1 


-Ja^Xm)n^Xm), 


m  =  l,2,...Nex+Nc. 
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(49) 

We  call  r"  the  displacement  residual  and  R™  the  traction  residual.  If  the  residuals  are  close 

enough  to  zero,  the  trial  set  of  u  and  p  represents  an  approximate  solution  of  the  boundary  value 
problem.  At  all  regular  boundary  points,  either  the  displacement  or  the  traction  vector  compo¬ 
nents  are  prescribed;  thus  only  one  of  Eqs.  (48)  and  (49)  would  be  used  in  assembling  the  overall 
system  of  equations.  For  points  that  are  in  the  cohesive  zone,  both  the  traction  vector  and  the 
displacement  discontinuity  vector  are  unknown  but  related  to  each  other  through  the  constitutive 
description  of  the  cohesive  zone  material  provided  in  Eqs.  (5)  -  (9).  From  Eq.  (49),  for  a  point  m 
that  lies  on  the  cohesive  zone,  the  residual  can  be  written  as  a  relation  between  the  displacement 
jump  vector  w(Xm)  and  the  traction  vector  p(Xm).  Note  that  points  on  the  crack  are  simply  spe¬ 
cial  cases  of  the  cohesive  zone  where  the  stiffness  has  decreased  to  zero  and  where  w(X"')  is  un¬ 
restricted  except  that  wn  (Xm )  >  0 .  Hence  the  complete  set  of  all  residuals  at  all  nodal  points  can 

be  calculated.  Note  that  in  order  to  incorporate  the  crack  opening,  crack  sliding,  and  crack  con¬ 
tact  modes  of  deformation  in  a  convenient  way,  field  quantities  in  the  discretized  equations  must 
be  transformed  into  the  local  orthogonal  coordinates  in  terms  of  the  normal  and  tangential  direc¬ 
tions.  We  now  turn  to  a  description  of  an  iterative  solution  scheme  for  solving  Eqs.  (48)  and  (49) 
for  u  and  p. 

For  a  node  m  on  the  external  boundary,  either  Eq.  (48)  or  Eq.  (49)  could  be  used,  depending 
on  whether  the  imposed  boundary  condition  at  that  point  is  a  displacement  or  a  traction  condi¬ 
tion.  For  the  iterative  solution  scheme,  in  general,  a  field  quantity  q™'l+]  at  the  (/+l)th  iteration  is 
written  in  terms  of  its  value  at  the  / 1,1  iteration,  q"J ,  and  an  increment  that  depends  on  the  re¬ 
siduals  as 

(50) 

where  R  is  the  appropriate  residual  for  the  field  quantity  q™  at  the  fh  iteration.  At  the  (/+l)th 
iteration  step,  the  displacement  component  at  this  node,  w”,/+1 ,  if  not  prescribed  as  a  boundary 
condition,  is  calculated  using  the  results  at  the  /th  iteration  by 

K'M  =  KJ  +  ®  ra  (u,p)/(cffff  +  h™),  (no  sum  on  a)  (51) 

where  co  is  an  adjustable  factor  of  relaxation.  The  displacement  residual,  ram  (u,  p),  is  computed 
using  the  nodal  values  at  the  {l+\),h  iteration  step  if  available,  or  at  the  Ith  iteration  step.  Simi¬ 
larly,  a  traction  component,  p™,M ,  if  not  prescribed  as  a  boundary  condition,  is  calculated  using 
the  results  at  the  Ith  iteration  by 

Pa,M  =  Pa1  +  fi>/?;(u,p)/(o.5-G;;),  (nosumon«)  (52) 

where  the  traction  residual,  R™  (u,  p),  is  dealt  with  in  the  same  way  as  ram  (u,  p).  For  a  node  m 

which  is  on  the  cohesive  zone,  we  calculate  the  displacement  discontinuity,  w",M ,  in  opening 
mode  by 

K  («,  p)/(*‘ - '  -  HZ ),  (no  sum  on  a)  (53) 
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where  m  is  another  adjustable  factor  of  relaxation,  and  tf"’1  is  calculated  using  Eq.  (31)  with 


w 


mj 


.  If  the  normal  component  of  w  at  the  (/+  \)th  iteration  step,  w",'+1 ,  is  negative,  it  is  set  to  be 


mj+1 


zero  for  no  penetration  of  the  crack  surfaces.  The  tangential  component  of  w  at  the  itera¬ 


tion  step,  <’,+’ ,  is  re-calculated  by 


<J"  -  ®  K(  u,  p)/(fJ  -  faP;-M  /  -H--  ),  (no  sum  on  r) 


(54) 


where  the  frictional  mode  of  the  cohesive  crack  is  activated,  w ^  if 


w™'1  <  wTd ;  otherwise 


= 

sary. 


.  Note  that  no  special  iterative  process  for  mixed  crack  opening  and  contact  is  neces- 


We  now  describe  the  marching  scheme  used  in  the  numerical  solution;  the  solution  procedure 
is  implemented  in  three  steps.  Suppose  that  the  solution  at  a  previous  loading  step  is  already 
known.  In  order  to  solve  the  problem  in  the  current  loading  step,  first  the  cracked  body  is  iterated 
to  be  in  equilibrium  with  all  the  current  cohesive  nodes  being  held  at  the  same  displacement,  us¬ 
ing  Eqs.  (51)  and  (52);  in  this  step,  the  cohesive  zone  material  law  is  not  used.  In  the  second  step, 
the  cohesive  zone  nodes  are  made  free  to  displace  according  to  the  constitutive  laws  described 
above.  The  whole  body  is  iterated  again  to  be  in  equilibrium  using  Eqs.  (51)  through  (54).  This 
holding  and  releasing  process  is  found  to  produce  quicker  convergence  for  the  iteration  process 
than  other  schemes  when  a  crack  is  advancing.  The  absolute  difference  of  either  displacement  or 
traction  at  a  boundary  node  between  two  next  iteration  steps  is  used  to  judge  the  convergence. 
The  best  relaxation  factors  of  co  and  m  can  be  obtained  through  some  trial  computations  and 
these  factors  vary  in  general  with  the  geometrical  configuration  of  the  body  and  loading  condi¬ 
tions.  Note  that  in  the  second  step,  the  fictitious  crack  tip  has  not  been  allowed  to  move.  After 
the  equilibrium  of  the  cracked  body  is  obtained  without  crack  advance,  in  the  third  step,  stresses 
at  a  point  of  the  fictitious  crack  tip  are  calculated  using  Eq.  (36);  note  that  a  constant  element  is 
always  used  at  the  fictitious  crack  tip  avoiding  the  difficulty  of  the  integration  due  to  the  hy¬ 
persingularity  in  the  fundamental  solution.  If  the  maximum  principal  stress  at  this  point  is  over 
the  critical  value,  i.e.  py,  the  fictitious  crack  tip  is  forced  to  advance  a  small  increment  in  the  di¬ 
rection  perpendicular  to  the  maximum  principal  stress.  This  is  accomplished  by  increasing  the 
length  of  the  last  cohesive  tip  element  if  the  advance  step  is  smaller  than  a  specific  value;  other¬ 
wise,  a  new  crack  element  would  be  added.  In  this  paper,  we  set  the  crack  advance  step  to  be 
equal  to  the  cohesive  zone  tip  element  size;  the  ratio  of  the  crack  advance  step  to  the  crack  ele¬ 
ment  size  could  play  a  role  in  determining  the  local  crack  patterns  in  situations  where  the  crack 
path  is  unstable,  but  this  issue  is  not  considered  in  this  paper.  Note  that  in  general,  a  numerical 
stepwise  scheme  of  crack  advance  may  generate  sharp  comers  between  the  approximating 
straight  elements  of  a  crack  which  may  concentrate  stresses;  these  stress  concentrations  are  of  a 
lower  order  than  the  crack  tip  concentration.  Investigating  the  role  of  these  comers  problem  is  of 
less  interest  than  that  of  examining  the  crack  tip;  thus,  in  the  present  paper,  we  neglect  the 
weaker  comer  singular  point  while  dealing  with  a  crack.  Note  also  that  at  the  current  loading 
level,  the  above  procedure  of  crack  advance  is  repeated,  and  the  iterative  solution  over  whole 
body  is  obtained  until  the  maximum  principal  stress  is  not  above  py  at  any  point  in  the  body.  This 
procedure  achieves  the  solution  in  the  current  loading  step.  The  load  is  then  incremented  and  the 
procedure  repeated  to  march  both  the  loading  and  the  crack  extension.  In  the  following,  we  shall 
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apply  the  method  formulated  above  to  some  problems  of  a  cohesive  crack,  and  demonstrate  the 
capacity  of  this  method. 

5.  Applications 

We  turn  to  a  demonstration  of  the  capability  of  the  boundary  element  strategy  for  handling  of  the 
cohesive  cracks  described  above.  In  order  to  explore  different  aspects  of  the  constitutive  model 
for  the  cohesive  zone  that  has  been  assumed  in  this  paper,  we  consider  two  problems.  The  first 
problem  concerns  a  single-edge-crack  in  a  rectangular  specimen  under  uniform  far  field  tensile 
loading.  For  small  initial  crack  lengths  the  crack  extension  in  this  configuration  is  unstable  under 
displacement  control;  however,  deep  cracks  exhibit  stable  crack  growth.  The  results  of  this 
problem  are  discussed  in  Section  5.1.  The  second  problem  considered  is  a  demonstration  of  the 
capability  of  the  formulation  to  handle  mixed-mode  crack  growth  and  crack  interaction;  using 
the  same  loading  as  in  the  first  problem,  but  with  an  offset  double-edge  crack  configuration,  the 
trajectories  of  the  two  edge  cracks  are  tracked.  The  results  described  in  Section  5.2  show  that  the 
model  duplicates  the  commonly  observed  behavior  of  two  approaching  cracks,  and  that  the  crack 
patterns  exhibit  a  sort  of  bifurcation  with  the  offset. 

The  following  conditions  describe  the  details  of  the  simulation  that  are  common  to  both 
problems.  The  domain  of  interest  is  a  rectangular  two-dimensional  region  of  length  l  and  height 
h,  taken  equal  to  /  in  these  simulations.  In  these  simulations,  all  length  quantities  are  normalized 
by  /;  while  this  is  not  a  natural  length  scale  for  the  fracture  problem,  it  is  convenient  and  easily 
interpreted.  If  it  is  desired,  a  scaling  to  an  intrinsic  length  scale,  such  as  the  critical  crack  opening 
displacement  for  the  cohesive  zone,  wy,  may  be  effected  easily.  The  deformation  of  the  specimen 
is  assumed  to  be  in  plane-strain.  An  initial  process  zone  of  length  0.02  is  assumed  to  exist  at  the 
tip  of  all  cracks  in  all  simulations  described  in  this  paper;  the  damage  parameter  Wd  is  assumed  to 
increase  linearly  from  zero  at  the  fictitious  crack  tip  to  vty  at  the  physical  crack  tip.  Note  that  this 
initial  damage  zone  can  be  prescribed  arbitrarily  and  is  in  general  unknown  in  the  physical 
problem.  The  boundaries  are  discretized  into  discontinuous  elements  with  a  constant  or  quadratic 
interpolation  over  the  element.  Specifically,  the  straight-line  segments  on  the  external  boundary 
are  divided  into  10  elements  with  quadratic  interpolation;  the  length  of  each  element  is  0.1  and 
the  nodal  spacing  is  0.0333.  The  initial  crack  (outside  of  the  cohesive  zone)  is  also  divided  into 
quadratic  elements  with  a  nodal  spacing  0.01.  The  initial  process  zone  is  divided  into  equally- 
spaced  elements  of  size  0.01  with  quadratic  interpolation,  except  for  the  element  at  the  fictitious 
crack  tip  which  is  taken  to  be  a  constant  element.  When  a  fictitious  crack  tip  advances,  one  new 
constant  element  of  size  0.01  is  added  as  the  corresponding  new  fictitious  crack  tip  element.  The 
old  fictitious  crack  tip  element  that  is  now  an  interior  element  within  the  cohesive  zone  is  then 
redefined  as  a  quadratic  element  using  three  nodes  inside  it.  The  constitutive  law  of  the  cohesive 
crack  is  assumed  to  be  represented  by  a  straight  line  as  shown  in  Fig.  10,  which  is  certainly  the 
simplest  one  defined  by  the  two  parameters  py  and  wf,  these  parameters  are  taken  to  be  py  =  0.01 
and  Wf  -  0.001.  Note  that  all  stress  and  traction  quantities  are  normalized  by  pi.  Poisson’s  ratio  is 
taken  to  be  0.3.  The  unloaded  specimen  is  initially  in  a  stress  free  state.  It  is  loaded  incrementally 
in  the  direction  perpendicular  to  the  top  and  bottom  boundaries  under  displacement  control.  The 
increment  of  the  loading  displacement  is  taken  to  be  IE-5  in  tension  and  -IE-4  in  compression. 
The  relaxation  factors  were  chosen  initially  by  trial  and  error  and  were  subsequently  fixed  at  co  = 
0.6  and  m  =  1.4. 
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In  the  following  simulations,  the  boundary  (including  the  cracks)  is  discretized  into  150  to  350 
nodes  depending  on  the  crack  length  resulting  in  a  system  of  300  to  700  equations  (or  degrees  of 
freedom).  In  addition  there  are  a  few  equations  that  provide  the  connection  between  tractions  and 
crack  opening  displacements  on  the  cohesive  zone.  If  the  stiffness  of  the  cohesive  zones  is  not 
modified  by  damage,  for  each  loading  step,  fewer  than  100  iteration  steps  are  needed  to  achieve  a 
solution  to  the  system  of  equation,  with  an  accuracy  of  IE-8  in  traction  components  and  IE-9  in 
displacement  components.  If  the  stiffness  of  the  cohesive  zones  changes  through  damage,  but 
still  without  propagation  of  the  fictitious  crack  tips,  slightly  more  than  100  iteration  steps  are  re¬ 
quired  for  each  loading  step.  If  a  new  cohesive  element  is  added  and  accumulates  damage,  each 
loading  step  requires  between  150  and  400  iteration  steps.  If  the  fictitious  crack  tip  movement  is 
sensitive  to  the  load,  quite  a  few  new  crack  elements  are  added  in  one  loading  step;  in  this  case, 
the  required  iteration  steps  may  be  as  large  as  1000,  depending  on  the  number  of  the  new  crack 
elements;  fortunately,  this  does  not  happen  often.  We  ran  the  code  for  the  following  problems  on 
a  Digital  Alpha  Workstation  (200  4/100).  Depending  on  the  size  of  the  discretized  system,  be¬ 
tween  10  and  30  iteration  steps  can  be  performed  in  one  second.  Run  times  are  as  short  as  Vi 
hours  for  small  crack  extensions;  if  long  crack  extensions  are  to  be  simulated  the  time  required 
increases  to  the  order  of  a  few  hours. 

5.1.  Single  edge  crack  under  tension 

A  number  of  simulations  were  performed  with  a  single-edge-crack  in  the  rectangular  geometry 
shown  in  Fig.  11.  The  crack  was  considered  to  be  along  the  line  of  symmetry  and  its  initial 
length,  ao,  was  varied  from  0.03  to  0.78.  The  variation  of  the  total  load  on  the  top  boundary  with 
the  imposed  elongation  obtained  from  the  simulations  is  shown  in  Fig.  12.  The  main  result,  of 
course,  is  that  in  these  simulations,  simply  by  prescribing  the  applied  loads  and  the  cohesive  law, 
initiation  and  growth  of  the  crack  appear  naturally.  These  results  are  the  mode-I  analog  of  results 
that  were  obtained  under  mode  III  by  Yang  and  Ravi-Chandar,  (1998)  using  a  finite  difference 
scheme.  A  number  of  remarks  regarding  these  simulations  are  listed  below: 

•  During  the  loading  process  of  a  stable  configuration,  the  fictitious  crack  tip  starts  to  advance 
first.  Note  that  the  critical  point  of  load  at  which  the  fictitious  crack  tip  starts  moving  is  of  no 
significance  since  it  is  dependent  on  the  initial  state  of  the  cohesive  zone.  With  continued 
loading,  the  physical  crack  tip  also  starts  to  move.  At  this  stage,  the  two  crack  tips  advance  at 
the  same  rate,  keeping  the  size  of  the  cohesive  zone-defined  as  the  distance  between  the  fic¬ 
titious  and  physical  crack  tips  -  a  constant  at  about  0.1.  This  is  considered  to  be  the  fully  de¬ 
veloped  equilibrium  cohesive  zone  size  under  this  condition;  of  course,  the  size  depends  on 
the  cohesive  material  model  and  the  geometrical  constraint  imposed  in  the  specimen. 

•  If  the  initial  crack  length,  ao,  is  small,  crack  initiation  is  unstable  even  under  displacement 
controlled  loading.  In  the  numerical  procedure,  this  is  manifested  by  the  fact  that  the  stress, 
calculated  at  a  point  ahead  of  the  fictitious  crack  tip  is  larger  than  py  even  after  the  crack  ex¬ 
tension  procedure  has  been  applied  over  a  length  of  many  cohesive  elements.  This  can  also 
be  inferred  from  the  fact  that  the  critical  load  for  the  onset  of  instability  decreases  with  in¬ 
creasing  initial  crack  length.  Thus,  crack  extension  and  structural  instability  coincide  in  these 
cases.  If  the  initial  crack  length  is  larger  than  about  0.48,  crack  extension  is  stable  under  dis¬ 
placement  control.  We  examine  later  some  issues  related  to  crack  paths  in  this  range  of  stable 
crack  growth.  The  behavior  indicated  by  the  sudden  drop  of  load,  for  initial  crack  lengths  less 
than  0.48  is  usually  referred  to  as  a  snap-back  instability  in  structural  mechanics.  Note  that 
the  snap-back  instability  is  also  predicted  by  the  stress  intensity  factor  based  linear  elastic 
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fracture  mechanics  and  is  simply  a  structural  feature  of  this  configuration;  it  is  not  influenced 
qualitatively  by  the  fracture  model. 

•  The  evolution  of  the  crack  opening  profiles  and  the  cohesive  tractions  at  the  onset  of  unstable 
crack  extension  are  shown  in  Figs.  9  and  10  for  a  few  initial  crack  lengths  in  the  unstable  re¬ 
gion.  Note  that  all  these  cracks  had  an  initial  cohesive  zone  of  length  0.02  prior  to  load  appli¬ 
cation.  For  long  initial  crack  lengths,  a  steady-state  cohesive  zone  of  length  0.1  develops. 
However,  for  small  initial  crack  lengths,  at  the  onset  of  unstable  crack  extension,  the  cohe¬ 
sive  zone  is  not  completely  developed.  In  other  words,  for  short  initial  crack  lengths,  struc¬ 
tural  instability  appears  before  development  of  a  self-similar  crack  tip  process  zone.  If  the 
dissipation  in  the  cohesive  zone  is  computed,  the  result  indicates  that  the  dissipation  is  a 
function  of  crack  length.  Note  that  the  data  used  for  Figs.  13  and  14  were  obtained  based  on 
the  numerical  procedure  of  the  loading  displacement  increment  equal  to  IE-5.  A  finer  load¬ 
ing  displacement  increment  may  produce  more  developed  cohesive  zones  representing  their 
critical  states  instead  of  the  ones  plotted  in  Figs.  12  and  14;  however,  this  characteristic  of  in¬ 
complete  development  of  the  critical  process  zone  at  the  onset  of  unstable  crack  initiation 
will  not  change.  A  similar  behavior  was  identified  in  the  mode  III  simulation  by  Yang  and 
Ravi-Chandar  (1998)  and  used  to  determine  the  range  of  applicability  of  a  single  parameter 
characterization  of  fracture. 

In  addition  to  the  structural  stability,  the  stability  of  the  crack  path  is  also  of  interest  in  crack 
problems.  The  single-edge-notched  geometry  used  here  is  very  stable  to  perturbations  in  the 
crack  path.  This  is  demonstrated  by  the  results  of  the  simulations  shown  in  Fig.  15.  Here  the 
crack  paths  obtained  from  three  simulations  are  shown.  The  only  difference  between  these 
simulations  is  that  the  edge  crack  is  offset  from  the  line  of  symmetry  by  d0ffset.  This  offset,  if  non¬ 
zero,  imposes  a  mixed-mode  loading  at  the  crack  tip  with  the  result  that  the  crack  path  is  no 
longer  straight.  The  path  stability  of  the  crack  is  indicated  by  the  fact  that  these  cracks  tend  to¬ 
wards  the  line  of  symmetry.  The  main  idea  here  is  to  demonstrate  the  capability  of  the  present 
boundary  element  formulation  to  track  arbitrary  crack  path  evolution.  Further  investigations  into 
the  influence  of  the  loading  on  the  crack  path  stability  are  still  under  progress. 

5.2.  Interaction  of  cracks  in  the  double  edge  cracked  specimen 

The  last  problem  we  consider  in  this  paper  is  a  double-edge-crack  loaded  in  tension  allowing 
crack  interaction.  This  problem  also  provides  an  opportunity  to  examine  cracks  that  exhibit 
structural  instability.  The  geometry  of  the  specimen  is  shown  in  Fig.  16.  The  initial  lengths  of  the 
two  edge-cracks,  a\0  and  ai0,  are  both  set  equal  to  0.28.  Each  crack  tips  is  again  provided  with  an 
initial  cohesive  zone  of  length  0.02.  In  order  to  examine  the  approach  of  the  two  cracks  towards 
each  other,  an  offset  d0ffset  is  provided  between  the  two  edge  cracks.  The  load-elongation  curves 
for  two  cases  of  d0ffset  equal  to  0.1  and  0.2  respectively  are  shown  in  Fig.  17.  The  deformed  con¬ 
figurations  showing  the  crack  paths  determined  from  the  simulations  are  presented  in  Figs.  18 
and  19.  The  main  observations  and  results  are  described  as  follows: 

•  Initiation  of  the  fictitious  and  physical  crack  tips  is  very  similar  to  the  cases  described  in 
Section  5.1.  Note  that  the  full  geometry  of  the  specimen  is  modeled  and  the  anti-symmetry  of 
the  offset  cracks  has  not  been  imposed  externally.  For  both  cases,  the  cohesive  ones  at  both 
crack  tips  are  fully  developed. 

•  For  both  simulations  with  d0gset  =  0.1  and  0.2,  after  initial  stable  extension  of  both  crack  tips 
over  a  length  of  about  0.1,  structural  instability  occurs  and  the  equilibrium  formulation  of  the 
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problem  may  no  longer  be  appropriate.  At  this  point,  the  cohesive  zone  is  fully  developed.  If 
we  assume  that  the  inertial  effects  are  small  in  dictating  the  crack  path,  we  may  be  able  to 
determine  the  path  evolution  of  these  interacting  cracks.  We  used  the  following  strategy  in 
our  simulations  to  determine  the  crack  extension  behavior.  At  the  onset  of  structural  instabil¬ 
ity,  the  trajectory  of  the  maximum  principal  stress  from  the  crack  tip  is  evaluated  first;  the 
cohesive  zone  -  i.e.,  the  fictitious  and  physical  crack  tips  -  is  then  extended  along  this  prede¬ 
termined  path,  at  fixed  global  displacement,  until  the  stress  ahead  of  the  fictitious  crack  tip 
falls  to  py  so  that  stability  is  restored.  After  propagating  the  crack  according  to  this  strategy 
for  a  length  of  about  0.1,  stability  is  restored.  We  note,  however,  that  this  scheme  did  not 
work  very  well  when  the  offset  was  very  small  indicating  that  inertia  effects  may  indeed  be 
significant  and  should  be  taken  into  account. 

•  It  is  interesting  to  note  that,  in  Fig.  18  for  d offset  =  0.1,  the  two  cracks  initially  ‘repel  each 
other,  pass  over  one  another  (overlap)  and  eventually  approach  each  other  along  a  curved 
path.  This  is  very  similar  to  the  experimental  observations  of  Melin  (1983).  On  the  other 
hand,  in  Fig.  19  for  d0ffset  =  0.2,  the  two  crack  tips  ‘attract’  each  other  right  from  initiation. 
From  a  different  point  of  view,  the  two  crack  tips  approach  the  line  of  symmetry  separately 
without  much  interaction  in  the  early  stages;  we  note  that  a  single  offset  crack  shows  the  ten¬ 
dency  to  approach  the  line  of  symmetry  as  demonstrated  in  Section  5.1.  After  initial  exten¬ 
sions  of  the  both  cracks  over  a  length  of  about  0.1  each,  the  left  crack  is  completely  shielded 
by  the  growth  of  the  right  side  crack  and  only  one  crack  grows  in  this  case.  The  change  of 
crack  patterns  with  d0ffset  is  a  bifurcation  phenomenon  of  significant  interest  (see  for  example, 
Mulhaus  et  al,  1996).  The  bifurcation  can  be  understood  in  the  following  terms:  consider 
unequal  extension  of  one  of  the  crack  tips  due  to  perturbations.  For  cracks  with  small  initial 
offset  distances,  perturbations  in  the  development  of  either  crack  tip  results  in  increased 
loading  of  the  lagging  crack  tip  while  for  large  initial  offset  distances,  these  perturbations  re¬ 
sult  in  an  increased  loading  at  the  leading  crack  tip. 

6.  Mixed-mode  fracture  with  a  cohesive  zone 

In  order  to  evaluate  the  influence  of  the  development  of  a  cohesive  zone  near  the  crack  tip  on  the 
mixed-mode  fracture  behavior,  the  formulation  described  above  was  used  to  simulate  the  ex¬ 
periments  described  in  Section  2.  The  geometry  is  as  shown  in  Fig.  11,  with  the  d0ffset  =  0  and  the 
load  applied  at  an  angle  #to  the  xj  axis.  The  results  are  described  in  Figs.  20  to  22.  In  Fig.  20, 
the  variation  of  the  load  against  the  load-point  displacement  is  shown.  For  comparison,  the  corre¬ 
sponding  linear  elastic  fracture  mechanics  results  are  also  shown.  Clearly,  the  two  results  indi¬ 
cate  a  good  comparison  suggesting  that  the  maximum  tangential  stress  criterion  should  be  a  good 
predictor  for  the  initiation  of  crack  extension.  Note  that  we  have  used  an  assumed  form  of  the 
cohesive  zone  model  and  hence  are  not  yet  able  to  make  quantitative  comparisons  with  the  ex¬ 
perimental  results.  The  comparison  of  the  crack  path  predicted  by  LEFM  and  the  cohesive  zone 
model  is  shown  in  Fig.  21  and  an  enlarged  view  near  the  initial  crack  is  shown  in  Figure  22.  The 
deviation  between  the  predictions  of  LEFM  and  the  cohesive  zone  model  in  determining  the 
crack  path  are  clear.  We  believe  that  this  is  an  indication  that  the  development  of  damage  near 
the  crack  tip  leads  to  significant  changes  in  the  evolution  of  the  crack  path.  Note  that  we  have  not 
provided  any  time  dependence  to  the  material  behavior  in  this  simulation;  this  will  introduce 
further  departure  in  the  development  of  the  crack  tip  process  zone  and  hence  in  the  crack  path 
evolution.  This  aspect  of  mixed-mode  crack  growth  is  under  further  study. 


19 


7.  Conclusions 


A  single-domain,  dual-boundary  integral  formulation  of  elastostatic  crack  problems  incorporat¬ 
ing  a  cohesive  zone  model  for  the  evolution  of  the  fracture  process  is  demonstrated  in  this  paper. 
The  cohesive  zone  is  modeled  as  a  damaging  material  with  a  prescribed  behavior  relating  the  ap¬ 
plied  force  to  the  crack  opening  displacement.  The  irreversible  nature  of  the  damage  is  intro¬ 
duced  through  the  damage  parameter  Wd,  the  maximum  displacement  experienced  by  a  point  on 
the  cohesive  zone  during  its  history.  The  introduction  of  the  cohesive  zone  necessitates  an  itera¬ 
tive  solution  procedure  to  solve  the  equations  resulting  from  the  boundary  integral  formulation; 
the  method  of  successive-over-relaxation  is  used  in  this  study.  In  terms  of  numerical  simulations, 
the  approach  described  here  presents  significant  advantages  over  grid-based  finite  element  meth¬ 
ods  since  the  present  formulation  (i)  does  not  force  development  of  the  crack  along  element 
boundaries  and  thereby  introduce  mesh  size  and  element  geometry  dependencies,  (ii)  does  not 
require  fine  nodal  spacing  except  within  the  cohesive  zone,  and  (iii)  does  not  require  remeshing 
with  crack  extension.  Three  example  problems  were  considered  to  demonstrate  the  capability  of 
this  formulation  to  handle  arbitrary  in-plane  crack  problems,  including  mixed-mode  problems, 
contact  problems,  and  crack  interaction  problems. 
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Figure  1.  Mixed-mode  loading  geometry  and  grip  attachments. 
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Figure  2.  Failure  interaction  curve  showing  critical  values  of  the  mode-I  and  mode  II  stress 

intensity  factors  at  crack  initiation.  The  prediction  of  the  MTS  theory  is  also  indicated. 
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Figure  3.  Dependence  of  the  mode-I  fracture  toughness  on  the  loading  rate. 


Figure  4.  Measured  variation  of  the  crack  kinking  angle  with  the  mode  mixity  paramter 
5  =  u/(1+li).  Predictions  of  the  MTS  theory  do  not  correspond  to  the  experimental  observations. 
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Figure  5.  (a)  Illustration  of  the  two  originally  coincident  points  of  a  cohesive  zone  in  the 

opening  mode,  connected  by  the  line  spring;  (b)  schematic  diagram  of  the  constitutive 
law  of  the  line  spring  in  terms  of  the  traction  and  displacement  jump.  Note  that 
instantaneous  loading  and  unloading  of  the  line  spring  is  given  by  the  slope  k. 
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Figure  6.  (a)  Illustration  of  two  originally  coincident  points  of  a  cohesive  zone  in  the 

contact/sliding  mode.  The  tangential  interaction  between  the  two  points  is  modeled 
using  a  Coulomb  frictional  force  and  the  cohesive  force,  if  the  cohesive  spring  is  not 
already  broken,  (b)  Sketch  of  the  variation  of  the  frictional  coefficient  of  the  contact 
surface. 
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Figure  9.  The  discretization  of  the  boundary  into  elements  is  shown  in  this  figure.  Each  element 
contains  one  or  mode  nodes  distributed  uniformly  within  the  element.  The  nodes  are 
internal  to  the  element,  indicating  discontinuous  elements. 


Figure  10.  The  two-parameter  constitutive  law  for  the  cohesive  zone  material  used  in  the 
simulations. 
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Figure  1 1 .  The  configuration  of  a  rectangular  specimen  with  a  single  edge  crack  of  initial  length 
cio,  under  displacement  controlled  elongation.  This  initial  crack  is  parallel  to  the 
loading  boundaries,  and  may  be  offset  from  the  symmetry  line  in  a  distance. 
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Figure  12.Load-elongation  diagram  of  a  set  of  simulations  with  a  crack  initially  lying  on  X2  =  0 
and  with  a  length  ao  from  0.03,  to  0.78.  An  initial  cohesive  zone  of  size  0.02  is  assumed  in  all  the 
simulations.  The  crack  growth  is  in  the  pure  mode  I.  For  ao  <  0.48  approximately,  crack 
extension  is  unstable  under  displacement  control. 
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Figure  13.  Crack  opening  profiles  at  the  critical  point  of  instability  in  two  of  the  simulations  for 
ao  =  0.03  and  0.18. 


Figure  14.  Traction  normal  component  p„  along  the  cohesive  zone  at  the  critical  point  of 
instability  in  four  of  the  simulations  for  ao  =  0.03,  0.08,  0.13,  and  0.18.  The  tangential 
component  is  zero  indicating  a  pure  opening  mode  for  the  cohesive  zone.  Note  that, 
for  very  short  initial  crack  lengths,  the  cohesive  zones  are  not  fully  developed. 
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Figure  15.  Crack  trajectories  in  the  simulations  with  different  offsets  of  initial  cracks.  Nodes 
along  the  newly-created  cracks  are  indicative  of  the  new  crack  element  ends. 


0 


original  crack 
newly  created  crack 


i — i — i — « — i — i — i — i — i — | — i — i — i — i — i — i — i — i — i — | — i — i — i — i — i — i — i — i — r 


Figure  16.  The  loading  configuration  of  a  rectangular  specimen  with  two  edge  cracks  on  the 
opposite  sides,  under  displacement  controlled  elongation.  The  cracks  are  initially 
parallel  to  the  loading  boundaries,  and  may  be  offset  from  the  symmetry  line  by  a 
offset  ^offset- 
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Figure  17.  Load-elongation  curves  in  the  simulations  of  crack  interaction  with  the  offset  d0ffset  as 
given  above.  Besides,  the  cracks  a\o  =  020  =0.28  initially,  and  each  crack  tip  was  assumed  to  have 
an  initial  cohesive  zone  of  size  0.02.  The  dashed  lined  indicate  the  loss  of  structural  stability. 
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Figure  18.  Deformed  configuration  of  the  specimen  with  two  edge  cracks  with  the  offset  dDffset  = 
0.1,  corresponding  to  an  extension  A  =0.0038. 
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Figure  19.  Deformed  configuration  of  the  specimen  with  two  edge  cracks  with  the  offset  d0ffset  = 
0.2,  corresponding  to  an  extension  A  =0.0038. 


Figure  20  Load-load  point  displacement  variation  for  mixed-mode  crack  growth,  showing  a 
comparison  of  the  predictions  of  LEFM  and  the  cohesive  zone  model. 
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Figure  21.  Comparison  of  the  crack  path  predicted  by  the  maximum  tangential  stress  criterion 
and  the  cohesive  zone  model. 


Figure  22.  Comparison  of  the  crack  path  predicted  by  LEFM  and  the  cohesive  zone  model;  only 
the  region  near  the  initial  crack  is  shown. 


